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1.  INTRODUCTION 


With  the  finite  element  method  recognized  as  an  efficient  tool  in 
solving  engineering  problems  in  solid  mechanics  during  the  1960's  the  applied 
mathematician  began  to  explore  the  finite  element  theory  in  search  for 
mathematical  foundations.  It  was  not  until  the  late  1960's  or  the  early 
1970's  that  the  finite  element  theory  was  rediscovered  as  a  part  of  the 
classical  approximation  theory,  functional  analysis,  and  Sobolev  space  theory, 
[  1 —5 ] .  In  recent  years,  finite  element  applications  into  fluid  mechanics, 
heat  transfer,  electromagnetic  field,  and  other  areas  have  been  increasingly 
active  [6,7]  as  the  concept  of  the  finite  element  theory  has  reached  a 
considerable  degree  of  maturity  at  least  in  linear  problems. 

In  fluid  mechanics  applications,  computational  difficulties  arise  from 
nonlinearity  associated  with  convective  terms  in  Navier-Stokes  equation  and 
energy  equation.  Complications  occur  also  in  the  potential  equation  governing 
high  Mach  number  flows.  The  standard  Galerkin  procedure  is  found  adequate 
for  low  Mach  numbers,  low  Reynold  numbers,  low  Peclet  numbers.  However, 
the  solution  stability  and  rates  of  convergence  in  the  Galerkin  formulation 
deteriorate  significantly  for  high  speed  flows  [7,8] .  Refinements  of  mesh 
would  improve  the  solution  [7,8  1  but  are  quite  costly.  A  recent  development 
in  an  effort  to  remedy  unstable  solutions  in  high  speed  flows  is  "upwind 
finite  elements"  [  9-12].  Unfortunately,  however,  upwinding  leads  to  un¬ 
desirable  damping  and  accuracy  is  not  guaranteed  [13,14],  To  overcome  these 
deficiencies,  more  recently,  the  optimal  control  least  squares  finite  elements 
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have  been  employed  with  all  second  order  derivatives  reduced  to  the  first 
order  resulting  in  a  symmetric  positive-definite  system  of  equations  [15,16, 
17].  The  penalty  method  is  utilized  in  handling  all  constraint  equations 
in  which  properly  adjusted  penalty  constants  play  a  major  role  in  acquiring 
convergence,  stability,  and  accuracy. 

The  main  part  of  the  present  report  consists  of  formulations  of 
Galerkin  finite  element  equations  applied  in  aerodynamics  with  emphasis  on 
transonic  flows.  In  particular,  we  develop  a  shock  element  with  continuous 
and  discontinuous  trial  functions.  Rankin-Hugoniot  conditions  are  ade¬ 
quately  satisfied  for  shock  elements.  Examples  are  shown  for  free  stream 
Mach  numbers  up  to  0.95.  Deterioration  of  solution  begins  to  appear  for 
M  >  0.95. 

Toward  the  end  of  the  contract  period,  the  author  proposed  a  new 
approach,  called  the  optimal  control  penalty  finite  elements,  in  order 
to  overcome  instability  in  solution.  Details  of  this  new  theory  are 
presented  and  numerous  example  problems  are  given  to  demonstrate  its 
validity  and  potentiality.  These  examples  are  compared  with  exact  solutions 
in  Tricomi  and  small  perturbation  equations.  Detailed  calculations  in 
the  full  potential  equation  using  this  new  theory  are  beyond  the  scope  of 
the  present  study. 
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2.  PRELIMINARIES  -  TRANSONICS 


The  main  difficulties  in  the  study  of  transonic  flows  stem  from  the 
nonlinearity  of  the  governing  equations,  existence  of  shocks,  and  the 
possibility  of  nonphysical  solutions  (entropy  condition  violated) .  These 
problems  arise  in  finite  difference  techniques  particularly  with  complicated 
geometries.  The  objective  of  the  present  study  is  to  introduce  the  finite 
element  analysis  in  an  effort  to  overcome  these  difficulties. 

We  consider  an  inviscid,  compressible  fluid  and  assume  that  the  flow 
is  potential  such  that 


V  •  pV  *  0 

p  -  00  f1 


in  ft 


1 

IV- 1 


(2.1) 

(2.2) 


V  -  V$  (2.3) 

where  <J>  is  the  velocity  potential,  p  is  the  density  of  the  fluid,  y  is  the 
ratio  of  specific  heats,  and  C*  is  the  critical  velocity. 

It  is  possible  to  rewrite  (2.2)  in  the  form 

P  -  (l  -  K  |7*|2)a  P0  (2.4) 

where  pfl  -  1,  K  -  ,  and  a  -  ~  . 

As  shown  in  Fig.  2.1,  the  flow  is  assumed  to  be  uniform  on  and 

tangential  at  T  ,  implying  that 

D 

-  V  •  n  on  T  (2.5) 

an  00 

-  0  on  fg 


3<J> 


(2.6) 
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In  two-dimensional  problems,  the  boundary  condition  for  the  trailing  edge 
may  be  chosen  as 

<J>  -  0  at  TE  (2.7) 

It  should  also  be  noted  that  is  discontinuous  across  the  cutting-line. 

The  lift  L  is  defined  as 

<J>+  -  <jf  -  L  (2.8) 

Here  we  must  use  the  Kutta-Joukowsky  condition  to  find  the  lift  such  that 

P+  -  P"  at  TE  (2.9) 

In  case  of  a  shock  the  Rankine-Hugoniot  conditions  must  be  satisfied 
such  that 

JpV  •  n]]+  -  JpV  •  nj  (2.10) 

and  the  tangential  component  of  the  velocity  is  continuous. 

Consider  a  disk  shown  in  Figs.  2.2  and  2. 3  .  Here  the  Kutta- 
Joukowsky  condition  is  automatically  satisfied  due  to  symmetry.  There 
are  two  possible  velocity  distributions  as  represented  by  Figs.  2.3  and 
2.4  .An  expansion  shock  prevails  in  Fig.  2.2  in  violation  of  thermodynamics 
laws  whereas  Fig.  2.3  is  physically  consistent.  It  is  necessary  to  impose 
the  entropy  condition  to  prevent  expansion  shocks.  It  is  well  known  that 
in  finite  difference  methods  the  entropy  condition  is  represented  by  an 
upwind  scheme  or  artificial  viscosities  in  the  supersonic  region.  In  the 
finite  element  method,  however,  the  Lagrange  multipliers  are  utilized  to 
impose  the  entropy  condition.  The  details  will  be  shown  in  the  following 


section. 
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3.  GALERKIN  AND  VARIATIONAL  METHODS 


3. 1  Galerkin  Approach 


When  viscosity  is  negligible  in  compressible  fluids,  the  analysis 
becomes  simplified.  Physically,  we  may  consider  the  fluid  to  be  a  perfect 
gas,  particularly  in  aerodynamics  design.  In  this  context,  we  begin  with 
Euler  and  continuity  equations  of  the  form 


V. 

J 


(3.1) 


epV 


(3.2) 


where  e  *  0  and  e  =  1  for  two-dimensional  flow  and  axially  symmetric  flow, 
respectively. 

A 

The  total  enthalpy  H  and  entropy  0  are  constant  along  each  stream¬ 
line.  Therefore,  we  must  have 
dH  =  H  V.dt  =*  0 

9  ^  *■ 

or 

H  V.  =  0  (3.3) 

9  i  1 


and 


dn  “  n’iVidt  =  ° 
or 

n,1v1  «  0  (3.4) 

Here  it  should  be  noted  that  (3.4)  is  valid  everywhere  except  for  passage 
through  a  shock  wave.  Along  the  streamline,  the  velocity  of  sound  is 
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9P 

3p, 


.fn 


n=const 

The  relation  between  pressure  and  density  in  an  ideal  gas  with  c 


constant  is  given  by 

„  P 


P  o 

—  ■  —  exp 
Y  Y 

P  Po 


ee)  ■ « -  ft) 


where  PQ,  p0,  and  n„  are  the  initial  conditions  and  y  =  c^/c^,  and 
function  C  remains  constant  along  each  streamline.  In  view  of  (3. 
(3.2),  we  obtain 

„  Pa2  Pa2 

P’i  =  3  P>i  -^’i  --f-  UnC)’i 

P 


Rewriting  (3.1)  in  the  form 

J  (V  V  ),  -  e...e.  V.V  +  i  P,.,  -  0 

2  j  j  i  ijk  kmn  j  n,m  p  i 

and  from  the  definitions  of  enthalpy  and  entropy, 


H  =  c  T 
P 


H  *  c  Jin  — 
V  PY 


1  P 


we  obtain 


Th , .  +  e,Jt  e.  v  v  -  (H  +  v  v  )  . 

i  ijk  kmn  jn,m  2  j  j  »i 


For  the  direction  normal  to  the  streamline,  (3.11)  becomes 


T  3n 


3H 

3n 


+  _ V,V  n 

ijk  kmn  j  n,m  i 


0 


where 

H  ■  H  +  i  V  V 
z  j  J 

For  two-dimensional  or  axially  symmetric  flow. 


(3.5) 

and  c 

v 

(3.6) 

the 

5)  and 

(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 


(3.12) 
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V  V  n  -  -  V(V  -  V  ) 
ijk  kmn  j  n,m  i  21  12 


(3.13) 


where 


V  =  V  n  -  V  n 
21  12 


Hence , 


(3.14) 


(3.15) 


V  -V  -  -  I  /ii  _  i n\ 

2,1  1,2  *  yan  yR  3n/ 

Multiplying  (3.1)  by  V  and  substituting  (3.5)  and  (3.2)  into  (3.1)  yields 

TiVi.i  +  rV-i-fl 

viVi,j-a>(vk  +  if)  -° 


(3.16) 


Expanding  (3.16)  gives 


fa  VI  r  V  V  /  \  eV 

1  -  l—  I  V  +  1  -  (-M  V - ^  (  V  +  V  )+  — i 

\a  /  J  1.  1  L  \a/J  2»2  a  \  »»2  2*1/  x2 


V  V 
1  2 


Adding  and  subtracting  — V4-  V 

a  l,  2 

yield 


(3.17) 

in  (3.17)  and  substituting  from  (3.15) 


2V  V 


eV 


- 7-2-  V  +  — -  *  E 

2  a  1,2  x 

2 


in  which 


E  -  /—  In  _  I  i£^ 

*  lyR  3n  a2  3ny 


(3.18) 


(3.19) 


This  is  the  most  general  form  of  the  governing  equation  for  two-dimensional 
or  axially  symmetric  inviscid  compressible  flows  for  subsonic,  transonic, 
and  supersonic  speeds. 


The  solution  of  (3.18)  for  a  two-dimensional  case  may  be  obtained  by 
the  method  of  characteristics  by  transforming  into  hodograph  plane.  This 
is  not  possible,  however,  for  axially  symmetric  flow  due  to  the  presence 
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I 

I 

I 

I 


\ 


I 

I 

! 

| 

I 


\ 

\ 


I 

! 

I 


I 

I 

I 

I 

I 

I 


of  the  term  eV  /x  .  If  the  entropy  r|  and  enthalpy  H  are  constants  In 

2  2 

the  direction  normal  to  the  surface,  then  the  quantity  E  in  (3.18)  vanishes. 
In  this  case,  (3.18)  written  in  terms  of  the  velocity  potential  takes  the 
form 


2  ii  92<E 

a7  3x  3y  3x3y 


(3.20) 


where 


a2  =  A  +  Bq2 

with  A  »  a2  =  a^  +  B<^,  subscripts  0  and  *>  indicating  the  stagnation 

point  and  undisturbed  state,  B  =  (1  -  y)/2,  and  q2  =  V2  +  V2. 

1  2. 

To  proceed  with  finite  element  formulations  from  (3.20),  it  is  con¬ 
venient  to  write  (3.20)  in  index  notation  as  follows: 

4>,ii  -  ~2  ”  0  (3.21) 

Setting  tp  =  we  obtain  the  local  Galerkin-f inite  element  equation 


'I2 


) .  <t> .  6 .  )  $  d£l 

i^  ij'  N 


*  0 


(3.22) 


Integrating  by  parts  with  first  derivatives  of  (f  held  constant,  we  get 


ANM^M 


G„  +  F„  +  Sv 


(.3.23) 


where 

*SH  ’  Ja* N.lVl*1 

S’  fai‘  WjVYj^mVq 

f  * 

S  ’  JT  ♦•i"iVr 

ss’-fTV 


C3 , 24a) 
(3,24b) 
(3.24c) 
(3.24d) 
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Here,  F^  represents  the  Neumann  boundary  conditions  on  solid  surfaces  or 
at  free  stream,  and  is  the  nonlinear  boundary  condition  on  the  surface  or 
pressure  discontinuity,  which  will  be  elaborated  in  Section  3.3.  For  a 
shockless  domain  we  set  »  0.  It  should  be  remembered  that  the  partial 
integration  of  (3.22)  for  the  nonlinear  term  is  not  exact  but  the  momentum 


is  conserved  more  exactly  than  the  case  in  which  this  term  is  not  integrated 


at  all,  commonly  known  as  nonconservation  form.  Note  that  two  of  the  non¬ 
linear  terms  can  be  integrated  exactly  by  setting 

and  (if)  “  -57  (y  (|f)  )  In  thls  case>  gn  and  SN  should  be  reduced 


J.  UUC  IIUU- 

(id)') 


by  2/3  of  the  integrated  values  corresponding  to  these  terms.  If  so,  the 


approximation  (<K^  held  constant)  affects  only  the  cross  derivative  term 
in  (3.20). 


For  Gn  *  *  0,  the  equation  is  reduced  to  an  ideal  incompressible 

flow.  The  assembled  form  of  (3.21)  can  be  solved  using  standard  techniques 
for  nonlinear  equations: 


theory  may  be  used: 

3V  3V  eV  2  V  3V 

(1  "  IT  +  !T+  —  “M«(1+Y)uI^r 
1  2  2  1 

where  U  is  the  free  stream  velocity  and  M  ■  U/a  is  the  free  stream  Mach 

CO 

number.  In  terms  of  potential  function,  we  may  also  write 


n  i  ^2<t>  i  e  -  m2  +  y)  34>  924> 

(1  ~  +  x^  M»  - U  ^  3# 


(3.26) 


This  equation  is  valid  for  transonic  flow.  If  the  right-hand  side  of 
(3.26)  vanishes,  then  it  is  valid  for  subsonic  and  supersonic  flows. 
The  pressure  coefficients  are  given  approximately  (higher  order 


terms  neglected)  by 


2V  eVz 

-V  -  -5*  °r  V 


\2 

2_  I  j 

■  U  3x  I  -  U2  \5x2/ 


(3.27a,b) 


The  finite  element  equation  corresponding  to  (3.26)  is  now  given  by 
the  matrices  for  two-dimensional  flow 


5  5=) « 

/(ti  "  'ir "  )  Vr 


(3.28a) 


(3.28b) 


(3.28c) 


(3 . 28d) 


where  dft  ■  dxdy.  Note  that,  in  contrast  to  (3.24),  and  are  exactly 


integrated  here.  The  solution  procedure  in  the  small  perturbation  theory 
is  the  same  as  in  (3.23),  but  with  a  considerably  simpler  nonlinear  term. 
The  jump  condition  across  the  shock  waves  is  satisfied  when  the  integrands 
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of  F^  and  are  equal  or 

I«  -  >0  £  -  *  <  fry’ll-  -  M -  ' 0 

where  o  denotes  the  jump  and  tij/n^  is  the  slope  of  the  discontinuity. 

In  general,  for  transonic  and  supersonic  flows  around  an  airfoil, 
shock  waves  are  expected  to  develop,  thus  leading  to  pressure  discontinui¬ 
ties  across  the  shock.  Because  of  the  presence  of  shock  waves,  solution  of 
the  equations  presented  here  will  suffer  a  breakdown.  A  remedy  of  this 
situation  as  well  as  general  discussions  of  shock  waves  will  be  given  in 
Section  3.3. 

3.2  Variational  Formulations 


The  variational  principles  for  compressible  flow  have  been  studied  by 
Bateman  [l8,19]  Herivel  £ 20 ],  and  Serrin  [2l] ,  among  others.  The  basic 
idea  of  deriving  the  variational  principle  in  nonlinear  equations  is  to 
determine  the  existence  of  symmetric  Frechet  derivative.  If  the  Frechet 
is  not  symmetric,  no  variational  principle  exists  in  general.  However,  if 
the  flow  is  irrotational  we  may  use  the  so-called  Bateman  principle  [l8,19] 
to  derive  the  variational  principle  for  compressible  flow. 


i  -  r  Pda  +  r 

*4  *4 


(3.29) 


where  P  *  A  +  BpT  with  A,  B,  and  y  being  the  constants.  The  first  variation 

of  the  integral  /  Pdft  attaining  the  maximum  is  equivalent  to  the  conserva- 
•'O 


tion  of  mass 


(pV'i 


0 


(3.29a) 
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To  prove  this  we  begin  with  the  relation  of  a2  »  3P/9p 

a2  -  B  Y  pY_1  =  a 2  +  (<£  -  q2) 

2a2 

y\2  00  2 

Denoting  q  ■  and  neglecting  the  constants,  the  variational  integral 

takes  the  form 


■  maximum 


The  first  variation  of  the  above  becomes 


With  some  algebra  it  can  easily  be  shown  that  the  first  variation  can  be 
simplified  to 


61 


i 


(q2  ~  1  4>»i54>,i  dQ  =  0 


or 

61 

Finally, 

61 


J ^  P4> ,  154> , ±  dn  *  Jcp604>, ±) , ±  -  6<J>(p4>,1),i] 


dfi  =  0 


J ^  P«5<1>  4>  *  ±ntdr  -  5<f)  (P4> ,  ±) ,  ± 


dfi  =  0 


Since  *  0  or  5<j>  *  0  on  the  boundaries,  we  now  have 

(P^)^  -  <PV*i  "  0 


Thus  the  validity  of  (3.29a)  h  been  confirmed. 

At  this  point,  we  wish  to  examine  the  relationship  between  the  variational 
principle  (3.29)  and  the  full  potential  equation  (3.21), 
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+  ~  G($)  ■  0  in  SI 


(3.30) 


with 


GW) 


t  [(«)' 


32*  ,  0  a<t>  34>  32<j>  ; 

3x  3x  3y  3x3y 


( iiV  a2$  1 

\3y J  3x3y J 


subject  to  the  boundary  conditions 
$  *  g(*,y)  on  T l 


i,ini  -  0  on  f2 


If  we  choose  y  *  2,  a2  *  q2  or  q2  =  a2 ,  and  keep  the  first  derivatives 

held  constant,  then  the  expression  (3.29)  can  be  integrated  by  parts 
as  follows: 


or 


61  (<t>)  «  ~  I2  ^j^j^’i^’i  *  ° 

<5I(<J>)  “  JTi^'±n±  ~  a2  ^’j^’i^’iHj^  ^  dF 

“  fj('i±  *  a2  ^i^j^ij)^  dn 


(3.31a) 


(3.31b) 


We  notice  that  the  integrand  of  the  last  term  of  (3.31b)  becomes  zero  as  it 
satisfies  (3.30).  Setting  (3.31a)  equal  to  (3.31b),  we  get 


5K<t>)  *  -  “2  dfi  +  J 


dr 


(3.31c) 


where 


f  *  ‘  ^ ’ ini  +  I2 


The  derivation  in  this  context  implies  that  the  Bateman  principle  (3.29)  does 
not  lead  to  the  full  potential  equation  (3.30)  unless  the  various  assumptions 
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as  proposed  here  are  valid.  Rewriting  (3.31c)  gives 

»  6  i  (*•!*•!  "  +j  f<*>  dFJ  (3. 31d) 

Thus  the  variational  principle  valid  only  under  the  assumptions  made  above 


takes  the  form 


I(«  ‘  j[(l  *fr 


(3. 31e) 


There  are  other  means  of  developing  the  variational  principle  in  an 
approximate  manner  using  Frechet  differentials  [22]  .  Toward 
this  end,  we  consider  the  nonlinear  problem 
N(u)  -  f  -  0 

and  the  variational  integral 


I(u,v) 


IvN(u)  -  vf  -  ug  dft 


The  first  variation  takes  the  form 


-  fj  6v  +  £n(u',v)  -  gjiu^ 


dfi  =  0 


This  may  be  called  an  adjoint  variational  principle.  An  alternative  approach 


is  to  write 


IJ 


N(u°)  -  f  5u  dfi 


in  which  variations  are  taken  with  respect  to  u  while  holding  u  constant. 
Upon  the  variation  we  set  u  ■  u°  and  recover  the  original  governing  equation 
as  the  "Euler  equation."  The  variational  principle  of  this  type  is  referred 
to  as  a  restricted  variational  principle  [22]  .  In  this  case, 
the  functional  is  not  stationary  unless 
N(u,u)  -  g  -  0 
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where  v  *  u  *  u°.  We  demonstrate  this  idea  for  the  compressible  flow 


below. 


The  derivation  of  a  restricted  variational  principle  for  (3.30)  begins 


6l(4>)  *  J  (4> * ii  -  “2  6<|>dn  =  0 


(3,32) 


With  the  first  derivatives  of  <}>  in  (3.32)  held  constant,  we  integrate  by 
parts: 

<5I((J>)  -  J  (}),ini6<J)dr  -  J  4>,i6(j),1dfi  -  l 

+  /  72  dfl 

Jn a  i  J  i  J 


<5I(<J>)  *  6  J  [4  (4*.^^  - 

-  JT  ♦•iVdr  +  Jr  P 


(3.33) 


The  variational  functional  to  be  made  stationary  for  the  solution  of  <p  assumes 


the  form 


I(4>)  *  fQ  i**i*’idfl"  j  I2 

-J  (fr.jn^dr  +  J  p  4,»i(>*j<>»inj^dr 


(3.34) 


The  integral 


l 


4>df  provides  the  Neumann  boundary  condition  on  solid 


surfaces  whereas  the  last  term  of  (3.34)  denotes  the  pressure  discontinuity. 
It  is  interesting  to  note  that  (3.34)  is  identical  to  (3.31e)  as  a  special 
case  of  Bateman  principle.  In  terms  of  the  finite  element  interpolation 
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functions,  we  may  write  (3.34)  in  the  form 

*  2  'WVm  ~  BNMPqVMq<(>N  "  FN^N  +  SN<*>N  (3.35) 

Minimizing  this  functional  with  respect  to  the  nodal  values  of  the  potential 
yields 

^Wm  "  BNMPQCt>M<M,Q  ‘  FN  +  SN  =  °  (3.36) 

Setting  the  second  term  equal  to  G^,  the  expression  (3.36)  is  seen  to  be 
identical  to  (3.23). 

The  finite  element  solution  of  the  flow  around  a  circular  cylinder 
at  low  Mach  numbers  formulated  from  the  pseudo-functional  of  the  type  (3.34) 
was  first  carried  out  by  Norrie  and  deVries  Q  23]  but  without  the  partial 
integration  of  the  nonlinear  term,  resulting  in  nonconservation  form. 

An  alternate  approach  proposed  by  Carey  f  24  ]  deals  with  the  stream¬ 
line  formulation  for  two-dimensional,  compressible,  inviscid,  irrotational 


flow.  The  governing  equation  is  given  by 

*’ii  +  M»(LT-i)]+  GW  ‘  ° 


(.3,37) 


C(40  -  ip 


’ii  W  ^  2  ^  +  ^’1  IF  ’i  =  ° 


(3.38) 
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Thus  the  variational  principle  (holding  constant) 

J  * 


1<W  'fa  H  f1  +  ^(^T1)]  ‘  J  ^  y-T1)  Ekj 
/r  j*'!”!*  -  U?(rT1)!kj£k£i>’i,>'jl!'-l,’l*j  dr 


'ki 


is  derived  as 

^’i^’i^’j^’J L  | 

(3.40) 


With  suitable  finite  element  interpolation  functions,  we  can  write  (3.40) 
in  the  form 


I(^)  2  Wm  *  BNMPqVmV^Q  "  Vn  +  Vn 


(3.41) 


The  corresponding  finite  element  equation  becomes 

Vk  -  WM  -  fn  +  ss  ■ 0  <3-42) 

with  the  appearance  of  this  equation  identical  to  that  for  the  velocity 
potential  formulation  (3.37). 

Carey  [  24  ]  ,  using  the  form  (3.37)  together  with  the  perturbation 
expansion  of  'p  instead  of  the  approach  (3.42),  presents  the  finite  element 
solution  of  the  fully-inf inite  compressible  flow  past  a  cylinder. 

Note  that  in  the  perturbation  solution,  the  Laplace  equation  is  solved 
first  and  the  further  solutions  repeated  with  effects  of  compressibility 
added,  the  basic  idea  being  similar  to  the  solution  of  (3.42)  in  which  the 
first  term  refers  to  the  incompressible  flow  whereas  the  second  te  repre¬ 
sents  a  compressibility  correction. 

3.3  Transonic  Aerodynamics 
3.3.1  General 

Transonic  gas  flows  are  governed  by  partial  differential  equations 
of  mixed  types.  The  theory  was  initiated  in  the  historical  memoir  of 


19 


Tricomi  ^25]  in  the  form 

which  is  hyperbolic  for  y  <  0  and  elliptic  for  y  >  0. 


(3.43) 


The  extreme  lateral  influence  produced  by  an  obstacle  in  a  stream  near 
Mach  number  unity  may  be  studied  by  the  linearized  equation  for  two-dimen¬ 
sional  flow  with  small  perturbations. 


(3.44) 


This  shows  that  when  =  1,  we  have  =  0,  and  therefore  the  perturbation 

velocity  depends  only  on  x.  This  implies  that  the  disturbance  pro¬ 

duced  by  the  body  is  propagated  laterally  with  undiminished  strength.  Since 
the  curvature  of  the  streamlines  remains  finite  as  y  increases,  it  follows 
that  infinite  pressure  differences  can  exist  between  the  surface  of  the  body 
and  points  at  great  lateral  distances  from  the  body.  Thus,  the  linear  theory 
invalidates  the  physical  conditions  and  we  may  resort  to  the  nonlinear  theory 


of  small  perturbation. 


*  0 

or  the  more  rigorous  and  general  equation  of  transonic  flow, 

/  2  3u  ,2  2%9v  _  3u  _ 

(a  ~  u  )  -g;  +  (a  -  v  )^  -  2uv  ^  =  0 

Written  alternatively, 

V2<p  -  G  =  0 


(3.45) 


(3.46) 


(3.47) 


c  i,f  M  a2*  +  2  i*  M  !!$_  + 1 mV  s2*  1 

a2  L V  3x  J  3x2  3x  3y  3x3y  y3y J  3y2  J 

For  unsteady  conditions,  the  small  perturbation  theory  provides 
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32ft  ,  32<p 

3x2  3y2 


2M‘ 


—  _  \ 
3t  \ 3x  / 


.2  92tf> 

00  3t2 


-  G 


(3.48) 


with 


m2  +  M2  U  +  y)  1$  9  jj> 

00  *■  00 


U  3x  3x2 

If  the  amplitude  of  airfoil  oscillation  is  small  in  comparison  with  the 
airfoil  thickness,  it  is  possible  to  linearize  the  unsteady  problem  by  taking 

A 

$  to  be  a  small  perturbation  on  <f>.  This  gives  the  following  linearized  dif¬ 
ferential  equation  for  <}>  in  nondimens ionalized  form: 


(1  -  Ml) 


2^  3  <f> 


15? 


-  M2  (y  +  1)  --  \  +  ifl  ,  e  3$  2M2— m2^  -  n 

»  ^  ’  9x  \3x  3x  / +  ip  +  ylj-  2M»lt  \lx /  0 


(3.49) 


with  x  =  x’/a,  y  =  y'/a,  t  =  t'U/a,  <J>  =  $'/Ua. 

For  the  flow  over  a  harmonically  oscillating  body,  we  set 


<J> 


I*1*] 


(3.50) 


where  $  is  a  complex  function  and  k  is  the  frequency  of  oscillation. 
Substituting  (3.50)  into  (3.49)  yields 

(1  -  m2>1!I  -  M2(Y  +  n  ii 1  +  lii  .  e  30  2,kM2  3$  , 

U  ^>95?  M»(Y  +  U  155  \l55  155/ +  Ip  +  7  37  "  2lkM»  155  + 

with  i  ■  /-T 


k2^^  =  0 


The  physics  represented  by  these  equations  is  predicted  reasonably  well 
qualitatively.  Some  significant  quantitative  results  from  wind  tunnel  ex¬ 
periments  are  also  available.  Shock  waves  are  an  essential  feature  of  tran¬ 
sonic  flow  and  their  appearance  on  a  moving  body  leads  to  a  rapid  increase 
in  drag  coefficient  with  increasing  Mach  numbers. 


I 
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The  analytical  solution  of  the  nonlinear-mixed  type  equations  is  a 
formidable  task.  In  recent  years,  significant  achievements  have  been  made 
by  workers  on  computational  fluid  dynamics  using  the  method  of  finite  dif¬ 
ference.  The  finite  difference  method  in  shock  wave  problems  was  studied 
by  von  Neumann  and  Richtmyer  [  26  among  others.  They  introduced  an  arti¬ 

ficial  dissipative  mechanism  so  that  the  shock  transition  would  become  smooth 
without  the  necessity  of  shock  fitting  [27]  in  which  the  incomplete  infor¬ 
mation  by  the  jump  conditions  must  be  supplemented  by  additional  infor¬ 
mation  coming  to  the  shock  from  the  fluid  behind  it.  Since  the  shock¬ 
fitting  technique  often  breaks  down  when  shocks  develop  spontaneously 
within  the  fluid,  von  Neumann  and  Richtmyer  avoided  the  direct  use  of  the 
Hugoniot  jump  conditions  but  retained  the  basic  conservation  laws  on  which 
the  Hugoniot  conditions  are  based,  in  which  the  jump  conditions  still  hold 
across  the  transition  layer.  Specifically,  the  dissipative  mechanisms  such 
as  viscosity  and  heat  conduction  are  included  so  that  the  surface  of  dis¬ 
continuity  is  replaced  by  a  thin  transition  layer  in  which  quantities  change 
rapidly,  not  discontinuously .  The  Hugoniot  relations  are  satisfied  such  that 
the  "smeared-out"  shock  travels  with  exactly  the  same  speed  as  a  discontinuous 
one  would. 

Alternative  forms  of  the  viscosity  term  have  been  proposed  by  a  number 
of  researchers  such  as  Lax-Wendrof f  [  28  } ,  MacCormack  [  29  ] ,  and  Murman 
and  Cole  [  30  ] ,  among  others.  In  the  recant  works  of  Jameson  [  31  ]  ,  the 
conservation  equations  are  centrally  differenced  for  subsonic  region  and 
upwind  differenced  for  the  supersonic  region  for  better  reprasentation  of 
shock  waves,  Wellford  and  Oden  [  32  ]  presented  the  finite  element  applications 
to  shock  waves  in  solids.  In  these  works,  the  derivative  of  the  axial 
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displacement  is  defined  in  terms  of  an  unknown  distance  at  which  a  shock 
front  travels.  The  finite  element  with  such  discontinuity  may  be  termed 
"moving  shock  element".  In  the  flow  field  around  a  body,  however,  a  shock 
is  at  least  two-dimensional.  Our  attention  is  focused  on  whether  there  is 
a  shock  discontinuity  at  any  given  location  and  time.  With  this  in  mind, 
we  introduce  here  a  special  finite  element  referred  to  as  the  "stationary 
shock  element"  p3]  . 

3.3.2  Element  Discontinuity  Method  for  Shock  Waves 

Let  us  assume  that  the  shock  wave  passes  through  an  element  causing  a 
pressure  jump.  Here  we  consider  an  explicit  discontinuity  at  the  center  of 
the  isoparametric  element  as  depicted  in  Fig.  3.1.  An  independent  inter¬ 
polation  of  <p  for  each  quadrant  is  given  by  (Fig.  3.2) 


- 

a !  +  a2  5 

+ 

ctp 

+  a4  Sn 

+ 

+ 

a6n2 

(3.52a) 

»(6)  - 

81+  82C 

+ 

e3n 

+  e4Sn 

+ 

+ 

*6n2 

(3.52b) 

*(Y)  - 

y  i + 

+ 

y  3n 

+  y„C  n 

+ 

y^ 

+ 

nn2 

(3.52c) 

6  j  +  62C 

+ 

s3n 

+  646n 

+ 

+ 

6sn2 

(3.52d) 

Here  we  introduce  a  total  of  24  constants  to  be  determined  uniquely.  Return 
to  the  equation  (3.52),  we  may  write  4  equations  for  <p  at  the  corner  nodes, 
8  equations  for  <f>  at  the  midside  nodes,  and  4  equations  for  <J>  at  the  center 
node,  resulting  in  16  equations.  We  obtain  8  additional  equations  for  the 
difference  between  the  slopes  of  4>  and  their  rates  of  change  at  the  center 
node.  They  are 


Figure  3.1  Isoparametric  Element 
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Figure  3.2  Representation  of  Discontinuous  Functions 
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3<{>(2)  34)(y)  „  3<})(y)  3({)(i5)  3<t>(<S)  34>(ot) 

2*  3 ~  '  W~  ,D3*TT  "  TT  ’  D“m  3n —  15n 


E 


3^(ct) 

TT2” 


32$(3) 

TH" 


„  _  32<J> 32(j)  Y  „  _  32<j^y*  32<(i^  „  _  32<}>^  32^a* 

2”  3n2  *"5tt2  ’  E  3m  Jjl  TT5-’  E|*m  “5n2  3ri2 

(3.53) 


With  these  24  equations  now  available,  we  solve  for  the  unknown  constants  to 
be  substituted  back  into  (3.52).  Thus  we  have 

<i>(o0  -  $<a)(5.n)<l>N  +  F^a)fjr,n)Dr  +  G^a)(5,n)Er 
<p(6)  -  ^6)(£,n)<|>N  +  F^s)(^,n)Dr  +  c£s)(s.n)Er 
D(y)  -  ^Y)(5.n)$N  +  F^Y)(5,n)Dr  +  G^Y)(^tn)Er 
<P(6)  -  ^5>(5.n)4>N  +  Fr6)U.n)Dr  +  Gr6)(5«Tl)Er 

where  N  -  1,2,..., 9  and  r  *  1,2, 3, 4.  It  is  clear  that  <t>^,  F^,  and  G^  represent 
the  interpolation  functions  for  continuity  of  and  discontinuities  of 
and  E^,  respectively  (See  Table  3.1). 

To  obtain  the  finite  element  analogue  of  (3.47),  we  assure  an  interpo¬ 
lation  field  for  the  velocity  potential  function  in  the  form 


(3.54) 


Here  V  represent  the  continuous  interpolation  functions  for  <f>  and  the  dis¬ 
continuous  interpolation  functions  for  derivatives  of  <J>  at  the  element  nodes; 

denote  the  nodal  values  of  <f>  plus  its  first  and  possibly  second  order 
discontinuous  derivatives.  An  orthogonal  projection  of  the  residuals  of 
(3.47)  onto  the  subspace  spanned  by  both  continuous  and  discontinuous  interpo¬ 
lation  fields  leads  to 


/  ($,,,  -  G)f  d£l  -  0 

Ja  11  m 


(3.55) 
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*7 
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*8 
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♦» 
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i+£n-S2-n2 

i-Cn-S2-n2 

*1 

i/2(S+S2) 

-1/2(C-C2) 

0 

0 

0 

i/2(m-n2) 

-i/2(n-n2) 

0 

'> 

0 

0 

1/2(C-C2) 

-l/2(5+S2) 

F- 

-i/2(n+n2) 

0 

0 

i/2(n-n2) 

Gj 

1/4(S+S2) 

1/4<S-S2) 

0 

0 

I 

g2 

0 

1/4 (n+n2) 

i/4(n-n2) 

0 

G3 

0 

0 

-i/4(C-£2) 

-i/4(S+S2) 

G4 

-i/A(n+n2) 

0 

0 

-i/4(n-n2) 

Table  3.1  Interpolation  Functions 
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which  is  identical  to  (3.22)  except  for  the  new  test  functions  ¥  replacing 

m 


the  standard  interpolation  functions  Integrating  by  parts  yields 


A  Q 
mn  n 


where 


G  +  F  +  S 
m  ni  m 


,dfl 


A  -  /  ¥  J  4d 
inn  l ^  m, i  n, i 

G  -  /  -2  ¥  .¥  .¥  ¥  ,dftQ  Q  Q 

m  a  n, x  p,j  q,i  m,j  TPp^q 

F  =  /  <f>,  ,n.¥  dr 

m  Jj,  i  l  m 

S  -  -  /  B.n  ¥  dr 
m  J  p  j  n  m 


(3.56) 

(3.57) 

(3.58) 

(3.59) 

(3.60) 


with  g  *  —2  cf>,  .<j>,  .<$>, . .  For  the  small  perturbation  theory,  these  matrices 
J  a  l  j  i 

are  of  the  form  (3.24),  now  equipped  with  shock  element  interpolation 


functions  instead  of  4. 


N’ 


The  boundary  conditions  are  given  by  F  represent- 

m 


ing  the  Neumann  type  on  V  and  denoting  the  surfaces  of  pressure  dis¬ 
continuity  on  T2  as  shown  in  Fig.  3.3.  We  introduce  a  notation  for  dis¬ 
continuity  given  by  (3.60): 

Sm  -  (S  )  -  -  (S  )  (3.61) 

m  m  1  m  c 

where  the  subscripts  1  and  2  refer  to  the  upstream  and  downstream  values 

of  S^,  respectively,  at  the  surface  of  pressure  discontinuity. 

From  the  element  geometry  (Fig.  3.1),  we  find  that  a  typical  17  x  17 

shock  element  influence  coefficient  matrix  A  takes  the  form 

mn 
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etc. ,  with  x  =  and  y  =  A^, 


A jj  being  the  standard  isoparametric  interpolation  function  derived  from 


x,y 


and 
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The  Gaussian  quadrature  integration  of  the  type  (3.62)  is  also  performed 
on  the  nonlinear  term  G^.  Thus  the  globally  assembled  finite  element  equa¬ 
tions  are  written  as 


Aa6Q8 


G  +  F  +  S 
a  a  a 


(3.64) 


Here  the  Neumann  boundary  conditions  F^  can  be  satisfied  automatically,  but 
we  require  special  treatment  for  the  shock  conditions  S^.  Toward  this  end, 
we  resort  to  direct  applications  of  Rankine-Hugoniot  conditions  (see  Fig. 
3.4)  to  obtain  an  equivalent  form  of  S^.  The  Rankine-Hugoniot  conditions 


are: 


piuni 
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together  with  momentum  normal  and  tangent  to  the  wave. 
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and  the  energy 

These  equations,  in 
equations  of  state, 
1  -  E  *  A 

with  A  ■  (1  -  £  j )  ( 1 


view  of  the  geometric  relations  and  the  appropriate 
lead  to  [34  J , 

-  d) 


Thus 


/u 


Y  ~  1 


d  -  tan’o/W’,  -1),  .J  (l  *  (y  -2l)(^)i’Z 

the  Rankine-Hugoniot  condition  results  in 

(l-A)u  -  u  “0 
ni  n2 


or 

[(1  -  A)  cos  a  -  t-a'nS^  l  -5)-  ]  »,  -  0  (3.65a) 

This  can  be  easily  recast  in  the  form 

q,  Q  “0  (3.65b) 

ka  a 

Here  u  t  is  the  resultant  local  velocity  upstream  of  the  oblique  shock.  The 

shock  angle  a  is  determined  from  the  discontinuity  values  at  the  element 

center.  Note  that  the  downstream  velocity  u2  is  eliminated  through  the 

deflection  angle  5,  which  is  determined  between  the  tangential  velocity 

u  and  the  downstream  velocity  u  .  The  shock  boundary  condition  matrix 
t2  *■ 

(3.65b)  is  now  equivalent  to  (3.61).  The  quantity  q^  is  called  the  shock 
boundary  matrix  with  k  denoting  the  number  of  shock  elements. 

To  enforce  the  shock  boundary  condition  (3.65b),  we  make  use  of 
Lagrange  multipliers  X  such  that 
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(3.66) 


To  replace  in  (3.64)  through  the  Lagrange  multipliers,  we  construct  an 


energy  functional  in  quadratic  form 
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+  F  .  Adding  (3.66)  to  (3.67),  we  obtain 
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Our  objective  is  to  seek  an  extremum  of  X  f°r  the  flow  field  satisfying  the 
jump  conditions  with  respect  to  every  and 

*  Iqh  5Qa  +  6  X  k  =  0 

3  X  _  „  ,  3  X 


For  all  arbitrary  values  of  6Qa  and  6)^,  we  must  have *  0  and  *  0, 
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Writing  (3.69) 

in  a 

compact 

form 

(3.69) 


BijXj  *  Yi 


(3.70) 


and  noting  that  Y  contains  R^  which  is  nonlinear,  we  may  cast  (3.70)  in  an 
iterative  form 

„(n) 


O  y(n+1> 
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Initially,  we  consider  a  shockless  domain,  with  Rq 

Vs  *° 


0: 


(3.71a) 


(3.71b) 


If  the  solution  (3.71)  yields  Q0  indicating  nonzero  D  .  then  (3.71)  will  be 

P  r 

replaced  by  the  expression  containing  the  jump  condition  (3.69)  in  the 
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subsequent  iterations.  In  this  iterative  scheme  R  will  be  kept  updated. 

The  magnitude  of  calculated  at  the  center  of  an  element  signifies 
the  strength  of  the  shock.  The  direction  of  the  shock  is  determined  by 
connecting  the  centers  of  shock  elements  (D^  4  0).  For  multiple  occurrences, 
interactions,  and  reflections  of  shocks,  it  is  advisable  to  start  with  the 
shock  element  influence  coefficients  applied  to  the  entire  domain.  At  each 
iterative  cycle,  we  remove  the  shock  boundary  matrix  when  is  found 
to  be  zero, 

A  qQq  -  R 
0t3  0  OL 

In  the  interest  of  simplified  iterative  calculations,  the  expression 


(3.55)  may  be  written  in  a  global  form  as 

fa  (Si  Vn  ' 0 


(3.72a) 


This  leads  to 


A  aQ.  =  G  +  F  +  S 
a0  0  a  a  a 


(3 , 72b) 


where 
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with  other  quantities  being  the  same  as  in  (3.64).  Notice  that  the  quantity 

K2  is  to  be  held  constant  during  each  iterative  cycle  and  updated  for  the 

following  cycles.  Had  we  used  the  Bateman  principle  with  y«  2,  then  we  have 

K2  «  (q/q)  with 
2a2 

a2  -  —  +  a2 

H  Y-l  » 
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Once  again,  the  shock  surface  boundary  condition  is  replaced  by  the 
Rankine-Hugoniot  condition  through  the  Lagrange  multipliers  (3.65).  Other¬ 
wise,  the  finite  element  equations  (3.72)  are  identical  to  those  derived 
from  (3.29)  or  (3.31). 

If  the  entropy  and  enthalpy  gradients  are  non-vanishing  normal  to  the 
streamlines,  then  we  need  to  modify  (3.55)  in  the  form 
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with  E  given  by  (3.19), 
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(3.74) 


where  the  angle  9  is  measured  between  the  x  axis  and  the  direction  normal 
to  the  streamlines.  In  this  case,  the  velocity  vector  is  given  by 
V  »  V<J>  +  7  *  <j> 

Here  4*  is  the  vector  potential.  This  situation  is  represented  by  E  ^  0  or 
the  entropy  and  enthalpy  gradient  vector- 
E  =  |  E'f  df  f  0 


E  =  /  E$  d 

m  Jt,  m 


which  is  added  to  (3.56): 


AQ-G+F+E+S 
mn  n  m  m  m  m 


(3.75) 


The  entropy  and  enthalpy  gradient  vector  is  calculated  from  the  relation 
between  the  pressure-density  and  the  local  velocity  distribution,  lagging 
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one  step  behind  the  iterative  solution.  It  is  observed  that  the  enthalpy 
gradients  in  general  are  negligible.  Direction  cosines  are  determined  from 
the  angles  between  the  x  axis  and  the  direction  normal  to  the  streamline 
which  corresponds  to  the  equipotential  line. 

The  17  x  17  element  coefficient  matrix  is  rather  large  in  size  ana  it 
may  be  advantageous  to  condense  it  to  a  smaller  size  via  standard  condensa¬ 
tion  schemes.  We  also  note  that  the  Gaussian  quadrature  integration  must 
be  performed  within  each  element  requiring  the  half  interval.  An  alternative 
approach  is  to  provide  an  independent  interpolation  field  for  each  quandrant 
so  that  the  Gaussian  integration  limits  run  -1  to  1.  In  this  scheme,  we 
provide  at  the  center  node  two  first  derivatives  of  <f>  for  each  quadrant, 
341^/35  and  etc.  Together  with  4  nodal  <j>  values,  there  are  6 

generalized  coordinates  for  each  quadrant.  Thus  the  interpolation  fields 
are  modified  to 
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Now  each  quadrant  can  be  treated  as  an  independent  element  having 
the  size  6x6.  In  the  assembly,  4>  is  common  to  all  quadrants  at  the  center 
node  whereas  all  the  second  derivative  terms  remain  as  generalized  coordi¬ 
nates,  independent  of  adjacent  quadrants.  In  this  approach,  although  the 
shock  strengths  are  not  directly  calculated,  the  definitions  given  in  (3.53) 
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Now  €:ach  quadrant  can  be  treated  as  an  independent  element  having 
the  size  6x6.  In  the  assembly,  <{>  is  common  to  all  quadrants  at  the  center 
node  whereas  all  the  second  derivative  terms  remain  as  generalized  coordi¬ 
nates,  independent  of  adjacent  quadrants.  In  this  approach,  although  the 
shock  strengths  are  not  directly  calculated,  the  definitions  given  in  (3.53) 
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enable  us  to  examine  the  existence  of  the  shocks  and  their  strengths. 

4 

3.3.3  Other  Approaches 

The  first  finite  element  applications  to  transonic  flows  are  those  by 
Chan,  Brashears,  and  Young  [  35  ]  who  employed  the  least  squares  in  conjunc¬ 
tion  with  the  Galerkin  approach  using  cubic  triangular  elements.  No  shock 
boundary  treatment  was  implemented.  Chan  and  Brashears  [  36  ]  further 
explored  their  method  to  include  an  unsteady  state.  Small  perturbation  theory 
is  utilized  in  their  studies. 

Subsequently,  Glowinski,  Periaux,  and  Pironneau  studied  the 
transonic  flow  by  optimal  control  [^7]  °f  distributed  parameter 
systems,  discretized  in  finite  elements.  The  method  of  steepest  descent 
was  used  to  solve  the  system  of  governing  finite  element  equations. 

Wellford  and  Hafez  [  38  ]  examined  a  mixed  variational  principle  for 
small  disturbance  transonic  flow.  The  interpolation  functions  for  the 
velocity  potential  and  the  velocity  were  introduced  in  the  minimization  of 
the  mixed  variational  principle,  resulting  in  two  sets  of  equations,  which 
may  be  solved  iteratively. 

Ecer  and  Akay  [  39  ]  studied  the  solution  of  full  potential  equation 
and  reported  difficulty  locating  a  shock  line  but  the  accuracy  of  the 
results  were  otherwise  good.  In  their  formulation,  the  Bateman  principle 
leading  to  (3.29)  was  used.  The  resulting  finite  element  equations  are 
similar  to  (3.72b)  but  no  shock  element  formulation  (3.52,  3.65,  3.69)  was 
considered  in  their  study. 


3.3.4  Unsteady  Transonic  Flow 


For  unsteady  conditions,  we  invoke  the  equation  of  the  form  (3.49)  and 
its  finite  element  analog, 

*Wm  +  SNM^M  "  “  FN  +  GN  +  SN  C3.77) 

with 

Sm  ■  fa  <Wa 
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Sm  =  /  2M2  <5  ft 

NM  00  N  dx 

For  oscillating  airfoil  problems,  we  return  to  (3.51)  with  an  interpolation 
field  of  the  form 


$  « 

<t  * 

1TN 

(3.78) 

(f)  * 

Vr)  +  %(!) 

(3.79) 

where  is  the  interpolation  function  and  and  are  the  real  and 

imaginary  parts  of  the  nodal  values  of  <p.  Substituting  (3.78)  into  (3.51), 
we  obtain  two  sets  of  equations  similar  to  (3.77)  but  one  for  real  part  and 
another  for  imaginary  part.  These  equations  will  be  solved  independently 
to  obtain  the  nodal  values  of  $  from  (3.79). 

In  general,  the  imposition  of  boundary  conditions  for  the  finite  ele¬ 
ment  solutions  is  simpler  than  in  the  finite  difference  method.  In  the 

0 

flight  conditions  of  an  airfoil,  the  exterior  boundary  is  infinite.  However, 
we  must  use  a  finite  geometry  along  which  the  far  field  analytical  solutions 
[40,  30]  may  be  matched  iteratively.  On  the  solid  boundary  of  a  body,  we 
require  that  the  velocity  normal  to  the  surface  be  zero;  i.e., 
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or 


34> 

3n 


b’ini 


Note  also  that  at  infinity,  the  perturbed  velocity  is  zero, 


V<{> 


The  imposition  of  Dirichlet  and  Neumann  boundary  conditions  follows 
the  routine  procedure  . 


For  unsteady  flows,  we  must  satisfy  additional  boundary  conditions, 
namely,  the  unsteady  Kutter  conditions;  that  is,  the  pressure  across  the 
trailing  vortex  sheet  must  be  equal. 


P 

where  C 
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(3.80) 

(3.81) 


Thus,  in  terms  of  potential  function,  the  unsteady  Kutter  conditions  are 

(3.82) 
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The  above  conditions  must  be  applied  along  the  upper  and  lower  surfaces  of 
the  mean  wake  position. 

If  the  body  is  oscillating,  then  we  require  that  velocity  on  the  sur¬ 
face  must  be  tangent.  Let  the  position  of  the  body  geometry  at  any  instant 
be  defined  as 

B(x,y,z,t)  =  0  (3.83) 


Then  the  surface  tangency  condition  can  be  satisfied  by  requiring  that  the 
material  derivative  of  the  function  B  with  respect  to  time  equal  to  zero 

[41]. 
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If  the  body  is  thin,  we  may  approximate  as 


3B  9I+1£3B+3£3B 

3t  3x  3y  9y  9z  3z 


(3.84) 


(3.85) 


For  a  two-dimensional  flow  (z*0) ,  if  we  let  B(x,y,t)  *  y  -  A(x,y,t), 


A<<  l,  ™  <  <  l. 


30  3A  3A 
■57  *  ^  +  9t 

where  the  function  A  is  taken  as 
A(x,y,t)  ■  ah(x)  +  bm(x,t) 


(3.86) 


C3.87) 


with  a  *  airfoil  thickness  ratio,  h(x)  =  thickness  distribution  along  the 
axis,  b  *  oscillation  amplitude,  and  m(x,t)  =  airfoil  oscillation  function. 


In  addition  to  these  special  boundary  conditions  (Kutter  condition  and 
oscillating  body  geometry),  we  must  satisfy  the  standard  boundary  conditions 
of  Dirichlet  and  Neumann  types  in  the  usual  manner. 


3.3.5  Example  Problems 


We  consider  here  a  missile  consisting  of  4  caliber  tangent  ogive  and 
9  caliber  afterbody,  which  was  studied  experimentally  in  a  wind  tunnel  by 
Spring  [  42  ]  .  (Fig.  3.5).  A  typical  finite  element  configuration  is  shown 
in  Fig.  3.6.  Note  that  the  finite  elements  tangent  to  the  body  are  installed 
with  two  sides  constructed  normal  to  the  body.  If  the  first  layer  is  made 
very  thin,  then  the  satisfaction  of  Neumann  boundary  conditions  will  be 
demonstrated  by  approxomately  equal  values  of  velocity  potential  between 
the  two  adjacent  nodes  in  the  direction  normal  to  the  body. 


Figure  3.5  Wind  Tunnel  Test  [32] 
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Figure  3.7a  Pressure  Coefficients  for  tl 
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Fig.  3.7c  Pressure  coefficients  for  the  problem  in 

Fig.  3.6,  circles  indicate  the  experimental 
results  [421,  M  »  0.95 
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The  solution  based  on  discontinuous  functions  (3.52)  together  with 
the  global  finite  elements  equations  (3.69)  is  obtained  through  the  iterative 
process  (3.71).  It  should  be  noted  that  less  computer  time  is  expended  if 
the  functions  (3.76)  are  used.  This  is  due  to  simpler  integration  ranging 
from  -1  to  1  instead  of  half-ranges  required  in  (3.62).  Both  integration 
schemes  appear  to  yield  identical  results. 

The  pressure  coefficients  for  M  =  0.85,  0.90,  0.95  are  plotted  in 
Figure  3.7.  Slight  deviations  between  the  experimental  results  and  finite 
element  calculations  are  assumed  to  be  due  to  the  far-field  boundary  condi¬ 
tions  of  finite  element  model  being  simulated  by  porous  wall  of  experiments, 
as  well  as  other  possible  factors  involved  in  both  computational  and  ex¬ 
perimental  procedures. 

Note  that  for  small  perturbation,  the  expression  will  assume 

a  simplified  form  (3.54)  with  all  other  procedures  remaining  the  same. 

3.3.6  Error  Estimates 

The  mathematical  error  estimates  on  mixed  equations  have  received 
growing  attention  in  recent  years.  For  simplicity  we  may  consider  the 
mixed  equation  of  the  type 

|j  IQfe.r>  £l  *  |p  ■  £  (3-s8) 

where  Q  changes  sign  in  ft.  Futher  simplification  leads  to  the  Tricomi  problem 

_  lljj  +  _  f  (3.39) 

y  Tf1 

Trangenstein  [43]  examines  Tricomi  equation  and  obtains  satisfactory  error 
estimates  on  the  elliptic  region  although  the  solution  on  the  parabolic  line 
y  •  0  and  in  hyperbolic  region  y<0  suffers  difficulty. 
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Subsequently  Aziz  and  Leventhal  [44]  studied  the  first  order  system 

Lu  »  Al|H  +  A2-|^-  +  A3u  ■  f  in  ft  (3.90) 

Bu  ■  0  on  T  (3.91) 

where  L  is  symnetric  positive  and  B  is  admissible  [45  ] .  Let  V  be  a  finite 

h 

dimensional  subspace  of  H1 (ft)  where  h  denotes  a  mesh  parameter  such  that 
and  4*^  £  V^.  Let  the  finite  element  equation  be  constructed  from 

(3.90)  and  (3.91)  to  find  U,  e  V,  such  that 


<UVVl 2<fl>  *  <NV*h>i!(a)  ■  (f-Vt ,(!!)•  V  \  *\ 

where  U  =  ♦  U.  and  ¥  ■  ¥  U..  This  leads  to  the  linear  system 

n  ii  h  li 


with 
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Aij  *  (L*i,fj)L2(a)  +  (NWi.2(n) 


(3.92) 


(3.93) 


F.  -  (f,^) 


It  can  easily  be  shown  [  46]  that 

liu  II  £  c  ||  f  II  c>0 

n  L2  (ft)  Lk  (ft) 

Now  let  the  subspaces  be  restricted  to  functions 
equal  to  k..  If  the  solution  u  in  (3.90)  and  (3.91) 
we  obtain  [46] 

||  u-U.  ||  -  0(hk)  (3.95) 

h  L2(fl) 

This  is  contrary  to  the  results  numerically  obtained  [44]  for  the  Tricomi 
problem.  A  similar  conclusion  arises  in  the  least  squares  approximation. 


(3.94) 

U.  of  degree  less  than  or 
h 

are  in  Hk+1 (ft)  O  C°(ft) ,  then 
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Consider  the  first  order  system 

0  <  x,y  <  1 

u(x,0)  =0,  0  <  x  <  1 

It  can  be  shown  [47]  that,  away  from  the  boundaries  x  *  0  and  x  ■  1,  the  error 
is  0(h2)  while  near  x  *  0  or  x  *  1  it  is  0(h).  The  overall  L2  error  is 

[//|u-uhr}“-0«,", 

one  half  less  than  optimal,  but  better  than  the  0(h)  predicted  by  LeSaint  [46] . 

The  theoretical  error  estimates  for  the  full  nonlinear  potential  equation 
are  not  available  at  this  time.  However,  our  numerical  results  indicate 
somewhere  between  0(h  )  and  0(h  )  pointwise  for  the  problems  formulated 
for  the  discontinuous  functions  described  earlier.  This  error  analysis  is 
for  the  potential  function;  the  rates  of  convergence  for  the  velocities 
calculated  from  the  potential  functions  are,  in  general  half  as  low. 

Let  us  now  investigate  the  rates  of  convergence  for  the  problem  of 
Fig.  3.6.  Upon  various  discretizations,  calculations  are  carried  out  to 
determine  the  slopes  of  the  pointwise  errors  at  x=10  on  the  airfoil  surface. 
The  results  for  M  =  0.9  are  presented  in  Fig.  3.8.  It  is  seen  that  the  rates 

of  convergence  for  the  full  potential  equation  and  small  perturbations  are 

.,.1.05,  .  ...1.13, 

C  (h  )  and  0  (h  ) .. 
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3.3.7  Remarks 

It  has  been  observed  that  the  Galerkin  Method  tends  to  be  unstable 
for  M  >  0.95,  although  stability  and  accuracy  may  be  restored  if  prohibitively 
large  number  of  elements  are  used.  Only,  recently,  the  author  has  investi¬ 
gated  the  optimal  control  penalty  least  square  finite  elements.  Following 
the  initial  success  of  this  method  in  convective  heat  transfer  problems  [17] 
we  began  exploration  into  the  Tricomi  equation  and  the  small  perturbation 
equation  [l5].  In  the  following  section  we  present  the  details  of  the  optimal 
control  penalty  least  square  finite  elements. 
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4.  OPTIMAL  CONTROL  PENALTY  LEAST-SQUARE  FINITE  ELEMENTS 


4.1  General 


[ 


We  consider  the  optimal  control  problems  of  the  type 
f  O20  +  G  ] 


Min 

4>e<!> 


(4.1) 


where  $  is  a  space  of  admissible  functions,  a  suitable  choice  for  $  being 
the  space  H*(fi).  An  alternative  form  of  (4.1)  may  be  constructed  by 
constraints 


94) 

u  =  9x 
94> 

V  =  -TT1- 

9y 


(4.2) 

(4.3) 


such  that 


Min 

4>e4>  *4  I 


r\B  + M- 


Here  X  and  y  are  the  penalty  constants.  It  is  seen  that  the  expression 
given  by  (4.4)  represents  the  first  order  system  of  integrodif ferential 
equation  and  leads  to  a  symmetric,  positive-definite  system  of  finite  element 
equations.  The  advantage  of  such  formulation  is  its  ability  to  provide 
type-independent  solutions. 

The  terns  of  G  for  Triconi  equation,  snail  perturbation  equation, 
and  full  potential  equation  are  expressed  as  follows: 

Tricomi: 

G  -  -  ,  a  -  1  -  y  (4.5) 


l 

r 
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Although  compact  in  form  as  seen  in  (4.10)  the  matrix  on  the  left 
hand  side  is  nonlinear.  The  solution  is  carried  out  using  the  Newton-Raphson 
procedure  together  with  conjugate  gradient  and  steepest  descent.  Let  us 
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consider  the  nonlinear  equations  resulting  from  (4.10)  in  the  form 


KijXj  '  Fi 


(4.11) 


where  i,j  denote  the  total  number  of  equations.  The  standard  Newton-Raphsor 
iterative  process  may  be  recast  in  the  form 


T(n)  (n+1)  T(n)  (n)  (n)  (n) 

ij  Xj  Jij  Xj  ‘  Ri  “  Bi 


(.12) 


where  is  defined  as 

o (n)  _  T(n)  v(n)  „(n)  v(n)  _(n) 
Bi  -  Jij  xj  -  Kij  xj  -  Fi 


(4.13) 


(4.14) 


It  is  convenient  to  calculate  the  Jacobian  for  each  element  and 
assemble  into  the  global  form.  The  solution  implied  in  (4.12)  may  begin 
with  the  conjugate  gradient  method  as  follows: 


V(D  .  F(0)  „  y(0)  =  _(0) 

Vi  Fi  "  KijXj  Fi 


a  =  ViVi  /  ViKijVj  =  Viri  / 
X<2>  -  x<l)+  a(1>V<l) 


(n+1)  ,  (n) 

i  i 


-  a(n)A  V(n) 

ij  j 


3  -  -  KijVj 

(n+1)  (n+1)  Q(n)  (n) 

Vi  ri  +  3  Vi 


(4.15) 


This  iterative  procedure  continues  until  the  residual  r^  10  .  The  values 

from  this  method  are  used  as  the  first  estimates  of  X^'s  in  (4.12), 
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In  order  to  solve  (4.12),  the  method  of  steepest  descent  is  used 
between  each  cycle  of  Newton-Raphson  iterations.  It  has  been  observed  that 
the  steepest  descent  is  more  efficient  than  the  conjugate  gradient  in 
stabilizing  and  inducing  a  rapid  convergence  for  the  Newton-Raphson  iterations. 
The  steepest  descent  equations  are 


(n+D  (n)  / _(n)  (n)  .  (n)  (n)  .  (n) 

\  \  +  Cri  ri  '  ri  Jij  rj)rk 


(4.16) 


where 

(n)  _  (n)  (n)  (n) 

ri  ~  Bi  "  Jij  Xj 


In  general,  approximately  10  cycles  of  Newton-Raphson  iterations  are  required 
for  satisfactory  convergence.  The  flow  diagram  shown  in  Table  4.1  outlines 
the  basic  steps  in  the  proposed  solution  procedure. 

The  treatment  of  shock  elements  can  be  handled  similarly  as  in 
Section  3.  However,  the  generalized  coordinates  in  (4.64)  are  replaced 
by  X^  and  the  jump  condition  (shock  boundary  condition)  equivalent  to  (4.66) 
is  written  as 


V*koXa  =  0 

These  observations  lead  to 

[Ka0  qka]  [XBl 

lqk6  0  J  Lvl 


(4.17) 


(4.18) 


Note  that  the  new  generalized  coordinates  Xg  contain  the  values  of 
derivatives  of  velocity  potential  but  not  the  differences  of  these  values 
between  two  adjacent  quadrants  of  an  element.  However,  the  shock  element 
can  be  defined  as  the  iterative  solution  converges. 
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4.3  Examples 
4.3.1  Trlcomi 


We  consider  the  solution  of  Tricomi  equation  in  the  form 
=  x4y  -  x2y4  +  fi 


(4.19) 


u  »  4x3y  -  2xy4 
v  =  x4  -  4x2y3  +  ^ 


which  satisfies 

vi!i  +  lii 

y  3x2  3y2 


0  <  x  <  1 


or 

■» 

We  let  the  penalty  constants  A  «  to  vary  between  small  and  large 
values  and  determine  the  optimum  magnitude.  As  seen  in  Fig.  4.1,  A  *  1.383 
appears  to  be  optimum.  Note  that  the  comparison  with  exact  solution  is 
quite  satisfactory  as  seen  in  Fig.  4.2.  The  rate  of  convergence  with  respect 
to  mesh  size  is  plotted  in  Fig.  4.3,  which  gives  approximately  0(h). 

For  y  <  0  the  Galerkin  method  or  the  standard  finite  difference  method 
fails  miserably.  The  present  approach  as  given  by  (4.1)  overcomes  the 
type-dependency  of  the  solution  and  achieves  an  excellent  stability  and 
accuracy  as  the  scheme  provides  a  type-independent  property  symmetric 
and  positive-definite  matrix.  The  accuracy  achieved  from  the  present 
technique  has  been  shown  to  be  far  superior  to  other  reported  results  [44] . 


Figure  4.1  Error  vs  Penalty  constants  for  mesh  size?  4x4, 

8x8,  12x12,  solution  of  Tricomi  equation.  Dirichlet 
Boundary  Conditions,  $  »  x4y  -  x2yu+  y  /21. 


»n||e|| 
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Hn  h 


Figure  4.2 


Error  vs  Mesh  size  at  optimal 
solution  of  Tricomi  equation, 


penalty  constant, 

<t>  -  x>  -  xV  +  y7/24- 


x  distance 


Figure  4.3 


Solution  of  Tricomi  equation,  comnarison  of  <f>  with 
exact  solution,  0  ■  x4y  -  xiy  +  y7/21. 


57 


4.3.2  Transonic  Equations 


Consider  the  cost  function  of  the  form 

J  *  /Ji ■  %  ' 0 " 5<*,y)) +  x  (u  *  H) 


+  v 


(4.20) 


where  G  is  given  by  (4.6)  or  (4.7).  For  small  perturbation  flow  we  define 


g(x,y)  »  -  (1-M*) (1-y) (2y-2x) (1-x)  +  2x(l-x)  + 

+  2M*  &&  [  -(l-y)xy  +  (l-x)(l-y)y]  [y(l-y)] 


(4.21) 


The  exact  solution  assumes  the  form 

<p  =*  (l-x)(l-y)xy  (4.22) 

The  optimum  penalty  constants  (A  *y)  can  be  seen  in  Fig.  4.4  and  the 
finite  element  solutions  are  shown  in  Fig.  4.5.  The  rates  of  convergence 
for  various  variables  are  plotted  in  Fig.  4.6.  The  stability  and  accuracy 
along  with  the  rate  of  convergence  achieved  in  this  method  are  superior 
to  any  previous  approaches  which  have  been  reported  to  date. 

With  conviction  that  the  optimal  control  penalty  least  squares  finite 
elements  are  capable  of  obtaining  type-independent  solutions  in  transonic  flows 
as  tested  with  exact  solutions,  the  procedure  for  shock  elements  and  dis¬ 
continuous  solution  can  be  combined. 

The  main  body  of  this  report  concerns  the  Galerkin  approach  as  pro¬ 
posed  in  the  original  contract  [48,49].  The  optimal  control  penalty  finite 
element  method  is  a  new  idea  which  has  grown  out  of  the  present  research. 

This  new  approach  warrants  extensive  further  research  as  the  author  is 
convinced  that  the  problems  of  discontinuity  and  exact  behavior  of  shock 
waves  can  be  most  effectively  studied  by  the  optimal  control  penalty  finite 


element  method. 


0 


0.5 


.0 


1 


1  .5 


Penalty  constant,  log1(y\ 

Fig.  4.4  Small  perturbation  transonic  eauation- 
Logarithmic  average  pointvise  error  vs 
penalty  constants.  $  »  (l-x)(l-y)xy 


I 


x  -  coordinate 


Small  perturbation  transonic  equation 
Comparison  with  exact  solution,  y-0.5 
8x8  mesh,  $  «  (l-x)(l-y)xy 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  work  carried  out  herein  represents  the  Galerkin  finite 
element  method  for  solving  aerodynamics  problems  with  emphasis  on 
transonic  flows.  A  shock  element  concept  was  proposed  and  computa¬ 
tions  have  been  carried  out.  In  this  process,  a  quadratic  iso¬ 
parametric  element  is  divided  into  quadrants  with  each  quadrant  having 
independent  trial  functions.  This  idea  allows  discontinuities  at  the 
center  of  an  element  and  shocks  are  allowed  to  develop  freely.  Rankin- 
Hugoniot  conditions  are  satisfied  accurately. 

Although  this  method  is  efficient  until  freestream  Mach  number 
reaches  approximately  0.95,  the  solution  seems  to  deteriorate  signi¬ 
ficantly  for  M  >  0.95.  Toward  the  end  of  this  contract  research 
period,  the  author  proposed  a  new  approach-optimal  control  penalty 
finite  elements.  This  method  is  suited  ideally  for  problems  of  dis¬ 
continuity  and  shock  waves  as  a  consequence  of  changes  in  the  type 
of  partial  differential  equations.  The  resulting  equations  are 
symmetric  and  positive-definite,  their  solution  being  type-independent. 
Numerous  examples  indicate  that  both  stability  and  accuracy  are  main¬ 
tained  very  satisfactorily  for  Tricomi  and  small  perturbation  equa¬ 
tions.  Detailed  calculations  applied  to  the  full  potential  equations 
using  this  approach  are  beyond  the  scope  of  the  present  study. 

In  summary  and  in  retrospect,  the  finite  element  method,  although 
extremely  successful  in  elliptic  boundary  value  problems,  has  just 
begun  in  aerodynamics.  In  due  time,  its  versatility  is  expected  to 
be  demonstrated  in  some  of  the  most  difficult  problems  in  aero¬ 
dynamics.  It  appears  to  this  author  that  a  full  exploration  of  the 
optimal  control  penalty  finite  elements  should  be  launched  since  this 
report  has  proved  its  potential. 
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